# This file produces figures in the text and supplementary information for
#   Goehring and Lowande's Public Responses to Unilateral Policymaking. 

# See the README for more information. 

# load required packages and data
require(tidyverse)
require(gridExtra)
require(kableExtra)
require(stargazer)

# plotting the simulation results
load('input-data/power-sim.Rda')

# round down for for viz
power$effect_size_t1[power$effect_size_t1==0.33]=0.32
power$effect_size_t1[power$effect_size_t1==0.36]=0.35
power$effect_size_t1[power$effect_size_t1==0.39]=0.38
power$effect_size_t2[power$effect_size_t2==0.17]=0.16
power$effect_size_t2[power$effect_size_t2==0.3]=0.29
keep=c(0.1,0.11,0.19,0.2,0.29)


f_power_lik_t1=ggplot(power[power$effect_size_t1%in%keep,]) + 
  geom_line(linewidth=1.2,aes(x=n,y=true.positive.likert.t1,color=as.factor(effect_size_t1))) + 
  scale_colour_manual(values=c('#d5ebff','#118bff','#00274C')) + 
  geom_hline(yintercept=0.7,color="#d86018") + 
  geom_hline(yintercept=0.8,color="#d86018") + 
  geom_hline(yintercept=0.9,color="#d86018") +
  scale_y_continuous(breaks=c(0.7,0.8,0.9),labels=c("70%", "80%", "90%")) +
  labs(x="N",y="Power",color='Effect Size') +
  theme_bw() + theme(legend.title=element_blank(),axis.title.y = element_blank())

f_power_dich_t1=ggplot(power[power$effect_size_t1%in%keep,]) + 
  geom_line(linewidth=1.2,aes(x=n,y=true.positive.dich.t1,color=as.factor(effect_size_t1))) + 
  scale_colour_manual(values=c('#d5ebff','#118bff','#00274C')) + 
  geom_hline(yintercept=0.7,color="#d86018") + 
  geom_hline(yintercept=0.8,color="#d86018") + 
  geom_hline(yintercept=0.9,color="#d86018") +
  scale_y_continuous(breaks=c(0.7,0.8,0.9),labels=c("70%", "80%", "90%")) +
  labs(x="N",y="Power") +
  theme_bw() + theme(legend.title=element_blank(),axis.title.y = element_blank())

f_power_lik_t2=ggplot(power[power$effect_size_t2%in%keep,]) + 
  geom_line(linewidth=1.2,aes(x=n,y=true.positive.likert.t2,color=as.factor(effect_size_t2))) + 
  scale_colour_manual(values=c('#d5ebff','#118bff','#00274C')) + 
  geom_hline(yintercept=0.7,color="#d86018") + 
  geom_hline(yintercept=0.8,color="#d86018") + 
  geom_hline(yintercept=0.9,color="#d86018") +
  scale_y_continuous(breaks=c(0.7,0.8,0.9),labels=c("70%", "80%", "90%")) +
  labs(x="N",y="Power") +
  theme_bw() + theme(legend.title=element_blank(),axis.title.y = element_blank())

f_power_dich_t2=ggplot(power[power$effect_size_t2%in%keep,])+ 
  geom_line(linewidth=1.2,aes(x=n,y=true.positive.dich.t2,color=as.factor(effect_size_t2))) + 
  scale_colour_manual(values=c('#d5ebff','#118bff','#00274C')) + 
  geom_hline(yintercept=0.7,color="#d86018") + 
  geom_hline(yintercept=0.8,color="#d86018") + 
  geom_hline(yintercept=0.9,color="#d86018") +
  scale_y_continuous(breaks=c(0.7,0.8,0.9),labels=c("70%", "80%", "90%")) +
  labs(x="N",y="Power") +
  theme_bw() + theme(legend.title=element_blank(),axis.title.y = element_blank())
  
  
g_legend=function(a.gplot){
    tmp=ggplot_gtable(ggplot_build(a.gplot))
    leg=which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
    legend=tmp$grobs[[leg]]
    return(legend)}
legend_1=g_legend(f_power_lik_t1)
legend_2=g_legend(f_power_lik_t2)


sim_plot <- grid.arrange(arrangeGrob(f_power_lik_t1 + theme(legend.position="none"), 
                         f_power_dich_t1 + theme(legend.position="none"), 
                         legend_1,nrow=1,widths=c(2.2,2.2,.6)),
             arrangeGrob(f_power_lik_t2 + theme(legend.position="none"), 
                         f_power_dich_t2 + theme(legend.position="none"), 
                         legend_2,nrow=1,widths=c(2.2,2.2,.6)),nrow=2)

ggsave('figures-tables/FIG-POWER-1.pdf', 
       sim_plot,
       width = 6,
       height = 5,
       units = 'in')



# Covariate balance ####



load('input-data/complete-data-v3.Rda')
covariates <- dt %>% 
  as_tibble() %>% 
  mutate(male = ifelse(sex == 'Male', 1, 0),
         female = ifelse(sex == "Female", 1, 0),
         high_school = ifelse(education == "High school or less", 1, 0),
         some_college = ifelse(education == "Some College or Vocational", 1, 0),
         college_grad = ifelse(education == "B.A. or B.S.", 1, 0),
         post_grad = ifelse(education == "Post-graduate or Higher", 1, 0),
         asian = ifelse(ethnicity == "Asian/Pacific Islander", 1, 0),
         black = ifelse(ethnicity == "Black", 1, 0),
         native_american = ifelse(ethnicity == "Native American", 1, 0),
         other = ifelse(ethnicity == "Other/Decline to State", 1, 0),
         white = ifelse(ethnicity == "White", 1, 0),
         twenty_five = ifelse(income == "1: Less than $25,000", 1, 0),
         fifty = ifelse(income == "2: $25,000 to $50,000", 1, 0),
         seventy_five = ifelse(income == "3: $50,001 to $75,000", 1, 0),
         hundred = ifelse(income == "4: $75,001 to $100,000", 1, 0),
         greater_hundred = ifelse(income == "5: More than $100,001", 1, 0),
         democrat = ifelse(partisanship == "Democrat", 1, 0),
         republican = ifelse(partisanship == "Republican", 1, 0),
         independent = ifelse(partisanship == "Independent", 1, 0),
         northeast = ifelse(region == "Northeast", 1, 0),
         midwest = ifelse(region == "Midwest", 1, 0),
         south = ifelse(region == "South", 1, 0),
         west = ifelse(region == "West", 1, 0))

# covariate balance in wave 1
cov_balance_w1 <- covariates %>% 
  dplyr::select(condition.w1, president.w1, age, male, female, high_school, some_college, college_grad,
                post_grad, asian, black, native_american, other, white, twenty_five, fifty, seventy_five, hundred,
                greater_hundred, democrat, republican, independent, northeast, midwest, south, west) %>% 
  group_by(condition.w1, president.w1) %>% 
  summarise(across(male:west, ~sum(., na.rm = T)),
            age = round(mean(age, na.rm = T), 1)) %>% 
  ungroup()

# transpose for printing
cov_balance_w1_t <- as_tibble(t(cov_balance_w1), 
                              rownames = "Variable",
                              .name_repair = 'unique') %>% 
  rename(V1 = ...1,
         V2 = ...2,
         V3 = ...3,
         V4 = ...4,
         V5 = ...5,
         V6 = ...6)

# turn sums into proportions
cov_balance_w1_t <- cov_balance_w1_t %>% 
  filter(!(Variable %in% c("condition.w1", "president.w1"))) %>% 
  mutate(across(V1:V6, as.numeric))

# pull out age since it is a mean not a count
mean_age_1 <- filter(cov_balance_w1_t, Variable == 'age')

cov_balance_w1_t <- cov_balance_w1_t %>% 
  filter(Variable != 'age') %>% 
  mutate(sum = V1 + V2 + V3 + V4 + V5 + V6) %>% 
  mutate(across(V1:V6, ~round(. / sum, 2))) %>% 
  dplyr::select(-sum)

cov_balance_w1_t <- rbind(cov_balance_w1_t, mean_age_1)

# make table
cov_balance_w1_t <- cov_balance_w1_t %>% 
  rename("Position, Obama" = "V1", "Position, Trump" = "V2",
         "Congress, Obama" = "V3", "Congress, Trump" = "V4",
         "Order, Obama" = "V5", "Order, Trump" = "V6") %>% 
  mutate(Variable = case_when(
    Variable == "male" ~ "Male", 
    Variable == "female" ~ "Female",
    Variable == "high_school" ~ "High school or less",
    Variable == "some_college" ~ "Some College or Vocational",
    Variable == "college_grad" ~ "B.A. or B.S.",
    Variable == "post_grad" ~ "Post-grad or Higher",
    Variable == "asian" ~ "Asian/Pacific Islander",
    Variable == "black" ~ "Black",
    Variable == "native_american" ~ "Native American",
    Variable == "other" ~ "Other/Decline to State",
    Variable == "white" ~ "White",
    Variable == "twenty_five" ~ "Less than $25,000",
    Variable == "fifty" ~ "$25,000 to $50,000",
    Variable == "seventy_five" ~ "$50,001 to $75,000",
    Variable == "hundred" ~ "$75,001 to $100,000",
    Variable == "greater_hundred" ~ "More than $100,001",
    Variable == "democrat" ~ "Democrat",
    Variable == "republican" ~ "Republican",
    Variable == "independent" ~ "Independent",
    Variable == "northeast" ~ "Northeast",
    Variable == "midwest" ~ "Midwest",
    Variable == "south" ~ "South",
    Variable == "west" ~ "West",
    Variable == "age" ~ "Mean age"
  )) %>% 
  rename("Covariates" = "Variable")

cov_balance_w1_t %>% 
  kbl(booktabs = T,
      format = 'latex',
      caption = "{\\bf Covariate balance in wave 1.} Displays the share of respondents in each treatment condition by demographic variable. \\label{tab:balance_1}") %>% 
  column_spec(2, width = "1.5cm") %>% 
  column_spec(3, width = "1.5cm") %>% 
  column_spec(4, width = "1.5cm") %>% 
  column_spec(5, width = "1.5cm") %>% 
  column_spec(6, width = "1.5cm") %>% 
  column_spec(7, width = "1.5cm") %>% 
  pack_rows("Sex", 1, 2) %>% 
  pack_rows("Education", 3, 6) %>% 
  pack_rows("Ethnicity", 7, 11) %>% 
  pack_rows("Income", 12, 16) %>% 
  pack_rows("Partisanship", 17, 19) %>% 
  pack_rows("Region", 20, 23) %>% kable_styling(font_size = 7) %>% 
  write_lines("figures-tables/cov-balance-wave-1.txt")

  
# Covariate balance in wave 2 ####
cov_balance_w2 <- covariates %>% 
  dplyr::select(condition.w2, president.w2, age, male, female, high_school, some_college, college_grad, post_grad, asian, black, native_american, other, white, twenty_five, fifty, seventy_five, hundred, greater_hundred, democrat, republican, independent, northeast, midwest, south, west) %>% 
  group_by(condition.w2, president.w2) %>% 
  summarise(across(male:west, ~sum(., na.rm = T)),
            age = round(mean(age, na.rm = T), 1)) %>% 
  ungroup() %>% 
  filter(!is.na(condition.w2)) 

# transpose for printing
cov_balance_w2_t <- as_tibble(t(cov_balance_w2), 
                              rownames = "Variable",
                              .name_repair = 'unique') %>% 
  rename(V1 = ...1,
         V2 = ...2,
         V3 = ...3,
         V4 = ...4,
         V5 = ...5,
         V6 = ...6)

# turn sums into proportions
cov_balance_w2_t <- cov_balance_w2_t %>% 
  filter(!(Variable %in% c("condition.w2", "president.w2"))) %>% 
  mutate(across(V1:V6, as.numeric))

# pull out age since it is a mean not a count
mean_age_2 <- filter(cov_balance_w2_t, Variable == 'age')

cov_balance_w2_t <- cov_balance_w2_t %>% 
  filter(Variable != 'age') %>% 
  mutate(sum = V1 + V2 + V3 + V4 + V5 + V6) %>% 
  mutate(across(V1:V6, ~round(. / sum, 2))) %>% 
  dplyr::select(-sum)

cov_balance_w2_t <- rbind(cov_balance_w2_t, mean_age_2)

cov_balance_w2_t <- cov_balance_w2_t %>% 
  filter(!(Variable %in% c("condition.w2", "president.w2"))) %>% 
  rename("Position, Obama" = "V1", "Position, Trump" = "V2",
         "Congress, Obama" = "V3", "Congress, Trump" = "V4",
         "Order, Obama" = "V5", "Order, Trump" = "V6") %>% 
  mutate(Variable = case_when(
    Variable == "male" ~ "Male", 
    Variable == "female" ~ "Female",
    Variable == "high_school" ~ "High school or less",
    Variable == "some_college" ~ "Some College or Vocational",
    Variable == "college_grad" ~ "B.A. or B.S.",
    Variable == "post_grad" ~ "Post-grad or Higher",
    Variable == "asian" ~ "Asian/Pacific Islander",
    Variable == "black" ~ "Black",
    Variable == "native_american" ~ "Native American",
    Variable == "other" ~ "Other/Decline to State",
    Variable == "white" ~ "White",
    Variable == "twenty_five" ~ "Less than $25,000",
    Variable == "fifty" ~ "$25,000 to $50,000",
    Variable == "seventy_five" ~ "$50,001 to $75,000",
    Variable == "hundred" ~ "$75,001 to $100,000",
    Variable == "greater_hundred" ~ "More than $100,001",
    Variable == "democrat" ~ "Democrat",
    Variable == "republican" ~ "Republican",
    Variable == "independent" ~ "Independent",
    Variable == "northeast" ~ "Northeast",
    Variable == "midwest" ~ "Midwest",
    Variable == "south" ~ "South",
    Variable == "west" ~ "West",
    Variable == "age" ~ "Mean age"
  )) %>% 
  rename("Covariates" = "Variable")

cov_balance_w2_t %>% 
  kbl(booktabs = T,
      format = 'latex',
      caption = "{\\bf Covariate balance in wave 2.} Displays the share of respondents in each treatment condition by demographic variable. \\label{tab:balance_2}") %>% 
  column_spec(2, width = "1.5cm") %>% 
  column_spec(3, width = "1.5cm") %>% 
  column_spec(4, width = "1.5cm") %>% 
  column_spec(5, width = "1.5cm") %>% 
  column_spec(6, width = "1.5cm") %>% 
  column_spec(7, width = "1.5cm") %>%
  pack_rows("Sex", 1, 2) %>% 
  pack_rows("Education", 3, 6) %>% 
  pack_rows("Ethnicity", 7, 11) %>% 
  pack_rows("Income", 12, 16) %>% 
  pack_rows("Partisanship", 17, 19) %>% 
  pack_rows("Region", 20, 23) %>% kable_styling(font_size = 7) %>% 
  write_lines("figures-tables/cov-balance-wave-2.txt")




# Correlates of Recontact ####











recontact_1 <- lm(usable.recontact ~ condition.w1 + president.w1 + topic.w1 + age + sex + education + ethnicity + income + partisanship + region, 
                  data = dt)

stargazer(recontact_1,
          out = 'figures-tables/correlates-recontact.txt',
          keep = c('condition.w1', 'president.w1', 'age', 'sex', 'education', 'ethnicity', 'income', 'partisanship', 'region'),
          omit.stat = c('ll','aic', 'ser', 'rsq', 'adj.rsq', 'f'),
          font.size = 'tiny', 
          omit.table.layout = 'n',
          digits = 3,
          star.cutoffs = NA,
          header = F,
          colnames = F, 
          dep.var.labels.include = F,
          model.numbers = F,
          column.labels = c("Recontacted"),
          covariate.labels = c("Congress treatment",
                               "Order treatment",
                               "Trump treatment",
                               "Age",
                               "Female",
                               "Some college or vocational training",
                               "B.A. or B.S.",
                               "Post-graduate or higher",
                               "Black",
                               "Native American",
                               "Other/Decline to state",
                               "White",
                               "\\$25,000 to \\$50,000",
                               "\\$50,001 to \\$75,000",
                               "\\$75,001 to \\$100,000",
                               "More than \\$100,001",
                               "Democrat",
                               "Republican",
                               "Midwest",
                               "South",
                               "West"),
          title='{\\bf Correlates of Recontact.} Reports OLS coefficients and conventional standard errors from models measuring the correlation between treatment conditions, demographic controls, and a binary dependent variable indicating whether the participant was successfully recontacted. The issue area treatment conditions are included in the model but excluded from the table (all topic coefficients have p-values greater than .05). ',
          label = 'tab:recontact')


# Attentive Respondents ####

attention_props <- tibble(check = c("First attention check", 
                                    "Second attention check", 
                                    "Both attention checks"),
                          prop = c(sum(dt$attent1 == 'Correct', na.rm = T) / sum(dt$usable.recontact),
                                   sum(dt$attent2 == 'Correct', na.rm = T) / sum(dt$usable.recontact),
                                   sum(dt$attentive, na.rm = T) / sum(dt$usable.recontact))
)

attention_props <- attention_props %>% 
  mutate(prop = round(prop, 2)) %>% 
  rename(`Attention check` = check, 
         Proportion = prop)

attention_props %>% 
  kbl(booktabs = T,
      format = 'latex',
      caption = "{\\bf Attentive Respondents.} Displays the share of recontacted respondents who passed the attention checks in wave 2. \\label{tab:attent_props}") %>% 
  row_spec(2, hline_after = T) %>% 
  kable_styling(font_size = 7) %>% 
  write_lines("figures-tables/attentive-respondents.txt")





# The correlates of attentiveness ####


attentive <- lm(attentive ~ age + sex + education + ethnicity + income + partisanship + region,
                data = dt)

stargazer(attentive,
          out = "figures-tables/attentive-correlates.txt",
          omit.stat = c('ll','aic', 'ser', 'rsq', 'adj.rsq', 'f'),
          font.size = 'tiny', 
          omit.table.layout = 'n',
          digits = 3,
          star.cutoffs = NA,
          header = F,
          colnames = F, 
          dep.var.labels.include = F,
          model.numbers = F,
          column.labels = "Attentive",
          covariate.labels = c("Age",
                               "Female",
                               "Some college or vocational training",
                               "B.A. or B.S.",
                               "Post-graduate or higher",
                               "Black",
                               "Native American",
                               "Other/Decline to state",
                               "White",
                               "\\$25,000 to \\$50,000",
                               "\\$50,001 to \\$75,000",
                               "\\$75,001 to \\$100,000",
                               "More than \\$100,001",
                               "Democrat",
                               "Republican",
                               "Midwest",
                               "South",
                               "West"),
          title='{\\bf Correlates of Attentiveness} Reports OLS coefficients and conventional standard errors from models measuring the correlation between treatment conditions, demographic controls, and a binary dependent variable indicating whether the participant passed attention checks in both waves.',
          label = 'tab:attentive')


